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Recent molecular dynamics (MD) simulations of liquid silica, using the "BKS" model [Van Beest, 
Kramer and van Santen, Phys. Rev. Lett. 64, 1955 (1990)], have demonstrated that the liq- 
uid undergoes a dynamical crossover from super-Arrhenius, or "fragile" behavior, to Arrhenius, or 
"strong" behavior, as temperature T is decreased. From extensive MD simulations, we show that 
this fragile-to-strong crossover (FSC) can be connected to changes in the properties of the potential 
energy landscape, or surface (PES), of the liquid. To achieve this, we use thermodynamic integra- 
tion to evaluate the absolute free energy of the liquid over a wide range of density and T. We use 
this free energy data, along with the concept of "inherent structures" of the PES, to evaluate the 
absolute configurational entropy Sc of the liquid. We find that the temperature dependence of the 
diffusion coefficient and of Sc are consistent with the prediction of Adam and Gibbs, including in 
the region where we observe the FSC to occur. We find that the FSC is related to a change in the 
properties of the PES explored by the liquid, specifically an inflection in the T dependence of the 
average inherent structure energy. In addition, we find that the high T behavior of Sc suggests that 
the liquid entropy might approach zero at finite T, behavior associated with the so-called Kauz- 
mann paradox. However, we find that the change in the PES that underlies the FSC is associated 
with a change in the T dependence of Sc that elucidates how the Kauzmann paradox is avoided 
in this system. Finally, we also explore the relation of the observed PES changes to the recently 
discussed possibility that BKS silica exhibits a liquid-liquid phase transition, a behavior that has 
been proposed to underlie the observed polyamorphism of amorphous solid silica. 

PACS numbers: 64.70.Pf, 66.20.+d, 65.40.Gr, 65.20+w 



I. INTRODUCTION 

Liquid silica is the archetypal "strong liquid," that is, a 
liquid whose viscosity ry and other measures of relaxation 
follow closely an Arrhenius behavior, log?/ ^ 1/T 
where T is the temperature. For most liquids, ?y in- 
creases significantly faster than an Arrhenius law as T 
approaches the glass transition temperature Tg-. these liq- 
uids are referred to as "fragile." 

Strong liquids such as silica are important as glass- 
forming systems. In a strong liquid, rj varies less rapidly 
with T near Tg, compared to a fragile liquid. As a conse- 
quence, a strong liquid can be held in a desired range of 
rj over a wider range of T than a fragile liquid. As every 
glassblower knows, this makes silica-based systems easier 
to manipulate just above Tg than any other commonly 
available liquid. 

The fundamental origins of strong behavior in glass- 
forming liquids is also a subject of continuing interest. 
We note in particular two recent developments. First, 
computer simulation work of Horbach and Kob P| using 
the "BKS" model of silica ^] has demonstrated that at 
high T, the model liquid exhibits fragile behavior, and 
then crosses over to a regime of strong behavior upon 
cooling. The work of La Nave and coworkers, based on 
instantaneous normal mode analysis, has shown that such 
a crossover is connected to a progressive reduction in the 



number of diffusive directions in phase space accessed 
by the system j^. Such a "fragile-to-strong crossover" 
(FSC) may be a general mechanism underlying the emer- 
gence of strong behavior, and has since been studied 
for a number of systems y\. We note that a crossover 
from a super-Arrhenius to Arrhenius dynamics may be 
a general feature of liquids around the so-called mode- 
coupling temperature Q, as is appearing to emerge in 
recent numerical studies, thanks to the the larger dynam- 
ical window made available by current computational 
power [3.l9lllCl||. However, the T region where equilibrium 
simulations can be performed is still limited, and does not 
allow for a precise statement of the T-dependence below 
the crossover temperature, as required to make final con- 
tact with models for the glass transition 0, 0, ^| . 

Second, a growing body of computer simulation re- 
search has established the importance of the potential 
energy landscape or surface (PES) for understan ding the 
dynamics of liquids near Tg laiMCllllllllIllIl 
llfl m m m M El. The pes refers specifically to 
U, the instantaneous potential energy hypersurface of 
the system, expressed as a function of the 3iV coordi- 
nates qt that specify the positions of the N atoms of the 
system; i.e. U = U{qi,q2, ■■■,q'iN)- The properties and 
topology of the PES have been carefully studied in the 
above cited works, predominantly in the case of fragile 
liquids, resulting in important i nsights into the equilib- 
rium 26] and out-of-equilibrium |27ll28l | thermodynamics 
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of supercooled states, and the connection between ther- 
modynamics and transport propertie |2lL |29|. However, 
the relationship of the PES to the dynamic properties of 
strong liquids is less well understood. In this paper, our 
focus is to clarify this relationship, and in particular, to 
determine if the FSC proposed for liquid silica can be 
connected to properties of the PES. 

Following previous studies of fragile liquids, our ap- 
proach is to apply the "inherent structure" formalism of 
Stillinger and Weber |l5| to molecular dynamics (MD) 
computer simulation data obtained for the BKS model of 
liquid silica. In this approach, the PES is partitioned into 
basins associated with the local minima of [Illlllll. 
Each minimum corresponds to a particular configuration 
of atoms and is called an inherent structure (IS). We 
denote by ejs the average potential energy of the IS's 
associated with the basins sampled by the equilibrium 
liquid at a given T and volume V. An IS and its energy 
can be obtained in computer simulation by carrying out 
a local minimization of U starting from an equilibrium 
liquid configuration. 

As we will describe in detail below, the evaluation of 
eis, combined with free energy calculations, allows us 
to calculate the configurational entropy, Sc, of the sys- 
tem Hill |23,|2l|2|. Sc determines the number of dis- 
tinct configurations explored by the system, in this case 
the basins of the PES. In a liquid, diffusion is associated 
with the exploration by the system of different basins of 
the PES. The work of Adam and Gibbs (AG) predicts a 
relationship (in the low T limit) between the character- 
istic relaxation time of the system and Sc The AG 
relation has been recently derived in a novel way |l2j . 
Generalizing the AG relation to the diffusion coefficient 
D, the AG relation can be written as. 
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where fiQ and A are presumed to be constant with respect 
to T. In the context of liquid silica, an interesting test of 
the robustness of the AG relation is possible by checking 
if Eq.n]is obeyed throughout the region in which the FSC 
occurs. If so, the AG relation then provides a basis for 
connecting transport behavior (quantified by D) to the 
properties of the PES (quantified by Sc and e/5). 

In a recent letter we showed that liquid BKS sil- 
ica behaves in a manner that allows Sc to be calculated 
from e/5, and that D and Sc are related as predicted 
by the AG relation. We were thereby able to show that 
the FSC in liquid silica is associated with a change in 
the T dependence of e/5, i.e. a change in the nature of 
the PES explored by the system as T decreases. We also 
found that this observation in turn has implications for 
other behavior observed in BKS silica, in particular, the 
possible occurrence of a liquid-liquid transition, and the 
behavior of the liquid as related to the so-called Kauz- 
mann paradox. 

To reach such conclusions, extensive MD simulations 



are required over a wide range of V and T, to calcu- 
late thermodynamic and transport properties, as well as 
careful examination of the IS properties. In addition, the 
absolute free energy of the liquid must be evaluated. In 
the present work, we provide a detailed description of 
the methods used to obtain the results summarized in 
Ref . |3]| , and also provide an expanded analysis and dis- 
cussion of the results. This work is organized as follows. 
In Section II we describe our MD simulations, includ- 
ing the interaction potential used. Section III provides a 
detailed description of the techniques we use to evaluate 
e/S: Scj and the total free energy of the liquid. Section IV 
presents the results of these calculations and provides a 
discussion of their implications. 



II. MOLECULAR DYNAMICS SIMULATIONS 

We carry out MD simulations at constant V. Most of 
our results are for a system of 444 Si atoms and 888 O 
atoms. A few simulations are carried out with a reduced 
number of particles (333 Si and 666 O atoms) in order to 
access longer physical times scales. Our MD simulation 
program is based on the "MDCSPC2" source code [s^ . 
We also reproduce a subset of our results using a code 
we have written independently of MDCSPC2. Note that 
all molar quantities are reported here in moles of atoms. 

Our model of atomic interactions in silica, denoted here 
as ^BKS, is based on the BKS potential, modified in two 
ways. First, the BKS potential energy for both the Si-0 
and 0-0 interactions diverges unphysically to negative 
infinity at sufficiently small distances, allowing "fusion" 
events to occur during simulation of high T systems. To 
prevent this, ^bks consists of the standard BKS poten- 
tial plus a short range term given by 
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(2) 



where r^j is the interatomic separation between an atom 
i of species /i, and an atom j of species v. To choose 
the parameters e^i, and cr^iy (see Tabled we first identify 
the value r.y — r*j at which the inflection of the stan- 
dard BKS potential occurs, below which the divergence 
to negative infinity develops. The parameters are chosen 
so that the new potential increases monotonically, and 
without inflections, as rij decreases for rij < r*j; and so 
that the difference between the new and the old poten- 
tials is small for ry > r*. . Similar approaches have been 
used in other works |33.l3l|. 

The second modification to the standard BKS poten- 
tial included in ^bks relates to the treatment of longer 
range interactions. As is common in implementations of 
the BKS potential, we calculate the long range contribu- 
tions to the Coulombic potential energy using the Ewald 
summation technique, with the dipole surface term set 
to zero |35|. The reciprocal space summation is carried 
out to a radius of six reciprocal lattice cell widths. In 
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this approach, the real space contributions to the BKS 
potential are usually cut off discontinuously at a specified 
distance, often chosen as L/2, where L is the length of an 
edge of the simulation cell. However, we study systems 
over a wide range of density p, and we desire a poten- 
tial for which the cut-off is independent of L. Also, for 
accurate determination of inherent structures, we wish 
to remove discontinuities in the potential energy arising 
from cut-offs, and to remove any L dependence from long 
range corrections associated with the cut-off. 

To achieve these goals, instead of discontinuously cut- 
ting off the real space potential contributions, we intro- 
duce a switching function. At a fixed distance = 
0.77476 nm the real space terms of the standard BKS 
potential are replaced by a 5th degree polynomial that 
tapers smoothly to zero over the range Rg < Vij < Rc, 
where Rc ^ 1 nm. The polynomial coefficients and the 
value of Rs are chosen so that the potential is continuous 



up to and including second derivatives at both 



Rf: 



and Tij — Rc] and so that the potential and its first two 
derivatives are monotonic for Rg < Vij < R^ and go to 
zero as Rc- These choices depend on the Ewald 

parameter a that occurs in both the real and reciprocal 
space contributions to the potential energy. For all L, we 
choose a = 2.5 nm~^ to ensure sufficient convergence of 
the potential energy in the reciprocal space summation 
for the densities studied. The value Rc = 1 nm where the 
switching function reaches zero is chosen to include third 
Si-Si neighbor interactions at most densities studied. 

The real space contribution to ^bks, denoted here as 
0, is therefore a piece-wise defined function of the form. 
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4>{Rs < nj < Rc) = Df^^{rij - Rs 



{nj > Rc) 



= 



Rs 
Rs 



(3) 



(4) 
(5) 



where erfc(x) is the complementary error function and e 
is the permittivity constant. The parameters are given 
in Table □ 

Note that the above modifications have the conse- 
quence that the average potential energy U and pressure 
P obtained using ^bks differ from those obtained using 
the standard BKS potential. We find that the differences 
are approximately independent of T along isochores. At 
T = 4000 K and p = 2.3072 g/cm^, we find that ^bks 
gives a P value 0.25 GPa greater than the standard BKS 
potential, and U is 2.5 kJ/mol higher. At the same T and 
p = 3.8995 g/cvp?, the respective differences are 0.9 GPa 
and 4.4 kJ/mol higher. These are not large differences 



on the scale of our measurements, and the qualitative be- 
havior of the system is, as shown below, consistent with 
that found in other studies based on the BKS model. 

For the free energy calculation to be described below, 
we also perform MD simulations using a binary Lennard- 
Jones (LJ) potential, in which two atomic species (also 
labeled "Si" and "O" ) occur in the same 1 : 2 proportion 
as in Si02. The LJ pair potential is of the form 
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The pair potential is cut off at rij 

is determined so that ^Lj{rij = 2.5s^^) = 0. These 
potential parameters are given in Tabled 

In order to obtain equilibrium properties we use the fol- 
lowing procedure. We equilibrate the liquid using veloc- 
ity rescaling for a time r long enough to allow Si atoms to 
diffuse an average of 0.2 nm, after significant relaxation 
of P and U have disappeared from the system history. 
The interval of velocity rescaling varies from 10 to 1000 
time steps depending on T. The time step for all runs 
is 1 fs, except for T = 7000 K, where the time step is 
0.5 fs. Velocity scaling is then turned off and the system 
is evolved in a constant {N, V, E) ensemble for at least 
IOt. {E is the total energy.) Using this approach, there 
is no appreciable drift in E during the constant (iV, V, E) 
phase, over which we calculate equilibrium quantities. 

For the lowest T where relaxation is slowest we modify 
this procedure to improve our sampling of phase space: 
we conduct up to five independent runs, with the con- 
stant {N, V, E) phase of each run lasting at least 2t. The 
reported properties (including T) are averages over both 
time and over the independent runs. Thus averages for 
low T state points are also calculated over a total of lOr, 
while at the same time the danger of an undetected trap- 
ping in an out-of-equilibrium state is reduced through 
comparison of the independent runs. 

The densities of the isochores simulated are given in 
Table ^1 while the state points studied are shown in 
Fig. n Note that we have studied the isochore at density 
2.3566 g/cm'^ in order to compare with previously pub- 
lished work ,3]. The simulations along this isochore are 
those that involve only 999 atoms; all others model 1332 
atoms. 



III. CONFIGURATIONAL ENTROPY 
CALCULATION 

In this section we calculate Sc from knowledge of e/5 
and the vibrational properties of the basins of the PES. 
Similar calculations have been carried out for water |2n| | , 
binary LJ mixtures ,21] and orthoterphenyl 0. 

We begin by writing the Helmholtz free energy F of 
the liquid along an isochore as [3 > 

F = eis{T) -TSc{eis{T)) + U^b{T,eis{T)). (7) 
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11- V 


Si-Si 


Si-O 


0-0 


standard BKS parameters 


A^^ (J xlO-'*^) 





28.845422 


2.2250768 


Bp,^ (nm"^) 





48.7318 


27.6 


C^^ (J nm« xlO"^^ ) 





-2.1395327 


-2.8038308 


^BKS' short range parameters 


e^. (J XIO-"") 





4.963460 


1.6839685 


a^v (nm) 





0.1313635 


0.1779239 


'I'flA's: switching function parameters 


(J/nm^ xlO-'^) 


-235.3529 


122.0161 


-53.16278 


E^y (J/nm'' xlO"''') 


-117.7993 


61.33742 


-26.25876 


F^. (J/nm^ xlO"") 


-23.83785 


12.33446 


-5.415203 


$Lj parameters 


(kJ/mol) 


23.0 


32.0 


23.0 


S;ii. (nm) 


0.33 


0.16 


0.28 



TABLE I: Potential parameters used in this work for both 
^BKS and Also required to specify ^^bks are: a = 

2.5 nm"\ Rs = 0.77476 nm, Rc ^ I nm, qsi = 2.4e and 
go = 1.2e, where e is the charge of an electron. 



Label 


V (cmVmol) 


V (cmVg) 


P (g/cm^) 


L (nm) 


A 


5.1359 


0.256443 


3.8995 


2.2479818 


B 


5.6423 


0.281722 


3.5496 


2.3195561 


C 


6.1486 


0.307012 


3.2572 


2.3869665 


D 


6.6550 


0.332292 


3.0094 


2.4507704 


E 


7.1614 


0.357577 


2.7966 


2.5114146 


F 


7.6677 


0.382863 


2.6119 


2.5692635 


G 


8.1741 


0.408147 


2.4501 


2.6246184 


H 


8.4984 


0.424340 


2.3566 


2.4157510 


I 


8.6804 


0.433426 


2.3072 


2.6777320 



TABLE II: Isochore volumes, densities, and simulation box 
sizes studied in this work. All runs model 1332 particles, 
except for isochore H, where we model 999 particles. The 
reference volume Vb corresponds to isochore I. 
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FIG. 1: State points simulated in the (a) V-D and (b) V-T 
planes. All points give results for equilibrated liquids. 



This expression separates F into two contributions: one 
stemming from the fact that the Hquid samples different 
basins of the PES, and one arising from the properties of 
the basins themselves. Sc is the entropy contribution re- 
sulting from basin degeneracy, i.e. it counts the number 
of basins associated with an inherent structure energy 
e/s 0, 01 • The vibrational part of the free energy f^b 
arises from the free energy of the basins. We note that 
the basin properties may change with e/s, e.g., the vibra- 
tional density of states of a basin associated with low e/s 
may differ from that of one with high e/5. In equilibrium, 
F is minimized with respect to e/5 and we obtain 



dF 



0=l-T 



dSc , dfyi 



IS 



deis deis 



(8) 



We will provide evidence below that the vibrational prop- 
erties of the basins do not change substantially from one 
IS to another. In this case dfvib/ deis = 0, and we may 
write, 



S,{T) = 5,(To) 



1 de 



IS 



To 



T' dT' 



dT'. 



(9) 



SciTo) is the configurational entropy at a reference T = 
To. 

Eq.O shows that the behavior of Sc{T) is controlled 
by eisiT). The work of Sastry and others E 113 
has shown that for fragile liquids, e/s decreases, and is 
negatively curved, as T decreases. In accordance with 
Eq. El Sc for fragile liquids also decreases, and is nega- 
tively curved, as T decreases. If liquid silica is fragile at 
high T, and then crosses over to strong behavior at lower 
T, we expect that on cooling, both e/5 and Sc will ini- 
tially behave as in a fragile liquid. However, from Eq. ^ 
strong behavior implies that Sc (and hence e/5) is con- 
stant with respect to T. Therefore, if the liquid passes 
from fragile to strong behavior as T decreases, the impli- 
cation is that eis(T) will decrease with T at high T, and 
then pass through a point of inflection, consistent with 
the approach to a constant at low T. 

To obtain Sc{Tq) we write. 



Sc{To) — S{To) — Svib{To), 



(10) 



where S{Tq) is the total liquid entropy and SyibiTa) is the 
entropy contribution arising from the vibrational proper- 
ties of basins in the PES. The vibrational part has both 
a harmonic and an anharmonic contribution, which we 
calculate separately: 



where l3f 



and 



Svib{T) — Sharva{T) + Sa7ih{T), 



^ 3N-3 / 
1=1 ^ 



(11) 



1^ h (12) 



Q — / 1 dEgnh 
iJanh — / rpf Qrj^i 0--I ■ 



(13) 



5 



Here Sharm is the entropy in the harmonic approxima- 
tion of the IS's obtained at a given T. The set {w^} 
describes the vibrational density of states of the IS's (de- 
tails below), h is Planck's constant over 2tt, R is the gas 
constant, and k is Boltzmann's constant. Eanh is given 
by 

Eanh{T) - E{T) - EharmiT) - eis{T), (14) 

where Eharm is the harmonic contribution to the energy, 
given by Eharm = 3i?T(iV - l)/iV. 

To obtain 675 we select 100 equilibrated liquid config- 
urations over the course of the MD run, perform a con- 
jugate gradient minimization '37] oilA^ and then average 
the results. As a stopping criterion for the conjugate 
gradient minimizations, we specify a relative tolerance of 
10~^ along line minimizations and a relative tolerance of 
10~^^ between line minimizations. In the case of isochore 
H our runs are the longest, and so we average over 1000 
configurations. 

In order to evaluate Sc and Sanh (using Eqs.ElandlH 
respectively) we first fit average values of e/s and Eanh 
to polynomials in T, and then evaluate the required in- 
tegrals analytically. The Eanh fit is constrained so that 
at T = 0, the value of Eanh and its first derivative are 
zero. This is consistent with Eanh being a correction to 
the harmonic approximation. 

It is important to recognize that the expression for 
Sanh in Ea. ll3l and hence the estimation of Sc via Eqs.|51- 
1141 is valid only under the assumption that the basin an- 
harmonicity does not change from basin to basin. To un- 
derstand this, consider the expression for Eanh in Eg. 1141 
from which Sanh is calculated. The terms contributing 
to Eanh are evaluated from equilibrium liquid proper- 
ties. Yet, an implication of Eq. E|is that an IS obtained 
from a liquid at (e.g.) 4000 K, when heated itself to 
4000 K, will give a value of Eanh equal to that obtained 
from equilibrium configurations at 4000 K. However, an 
IS obtained from a liquid at 3000 K, when heated to 
4000 K will not necessarily yield the value of Eanh found 
from equilibrium configurations at 4000 K, because IS's 
obtained from different T may be in basins of different 
shape, and hence different anharmonicity. If the basin 
shape does change with e/s, then Sc in Eq. |51will be in- 
fluenced by an additional contribution. Moreover, Eg. 1141 
would be invalid and Eanh would have to be obtained in a 
different way, possibly by a careful heating of individual 
basins obtained from the equilibrium liquid at different 
T. Such heating experiments must be performed with 
care, as the system must not diffuse out of the basin if 
accurate results are to be obtained. 

To test if basins associated with different T have differ- 
ent shapes, we carry out runs in which IS's from different 
T and V are rapidly heated. We find that Eanh is the 
same for all basins belonging to the same isochore up to 
high T (Fig.Hl). 

Based on the relations justified above, we can evaluate 
Eanh{T) and Sanh{T) from a knowledge of eis{T). We 
can also evaluate Sc(T), up to a constant, from eis{T). 
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FIG. 2; Test of T-dependence of basin shape. IS's from three 
T, and for three isochores, are rapidly heated in order to 
confine the sampling to a single basin. Using velocity scaling, 
T is increased from K to 7000 K over 100 fs. Each curve 
is an average over 10 runs. The curves for the same isochore 
are approximately the same, indicating that the anharmonic 
contributions to the vibrational energy can be assumed, for 
the present purposes, to be the same for each basin 



To complete an evaluation of Sc{T), we need to estimate 
both Sharm{T) and S{Tq) for each isochore to be studied, 
as described in the following two subsections. 



A. Harmonic Entropy of Inherent Structures 

We define Sharm of the liquid as the average harmonic 
entropy of IS's sampled from the liquid. When a liq- 
uid configuration is quenched to its corresponding IS, it 
becomes a mechanically stable solid, and is to a first ap- 
proximation, harmonic. To calculate the entropy of an IS 
in the harmonic approximation, we require its vibrational 
density of states. As each IS is an atomic configuration 
at a local minimum of the PES, we expand the expression 
for lA about the local minimum: 



U ^ eis + ^^qi 

1=1 j=i 



dqidqj 



(15) 



Here, the set {qi} specifies the 3N atomic coordinates, 
and the notation "g = 50" denotes that the second deriva- 
tives are evaluated at the minimum energy configuration. 
We then define a Hessian matrix. 



1 



(16) 



9=90 



where is the mass of the atom associated with coor- 
dinate qi. Since the system is at a minimum, Hij has 
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FIG. 3: f2 as a function of T for (a) isochore A, (b) isochore D, 
and (c) isochore I. Also shown is the standard deviation about 
the mean value based on 100 samples. We note that a differ- 
ence in Q. of 0.01 yields a change in entropy of 0.24 J/mol K. 
For the purpose of this work, we therefore consider to be 
constant along isochores. (d) Q.{T = 4000 K) as a function of 
V\ this is the value of used in our calculations. 
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FIG. 4: S7 as a function of eis for (a) isochore A, (b) isochore 
D, and (c) isochore I. The error bars represent the standard 
deviation about the mean value based on 100 samples. 



eigenvalues {hi\ all greater than zero, except for three 
zero eigenvalues which account for the three independent 
translations of the entire system. These three eigenvalues 
are excluded in calculating the harmonic entropy. The Wi 
appearing in Eq. 1121 are defined as Wi — y/hl- We note 
that a particular Hessian matrix corresponds to an IS 
obtained from a liquid configuration at a certain T. It is 
this T that we use in Eq. ^1 

We find, perhaps surprisingly, that the spectrum of Wi 
does not change appreciably with T along isochores. We 
plot in Fig. 121 the quantity. 



n{T) = 



^— y 

iN-'d ^ 

i—1 



(17) 



to show that the contribution to Sc from changes in Q. as 
T, and hence e/g, is varied is negligible. This quantity 
is part of the expression for Sharm(T). It captures the 
average quadratic shape of a basin and hence determines 
any dependence of Sharm on e/s. The plot shows fl not 
to vary appreciably with T, and we conclude that there 
is no contribution to Sc from /t,i&, at least not from the 
harmonic portion. 

To confirm this approximation, we show in Fig. 0] the 
variation of f2 with ejs- We find that the change in Q can 
be as large as Ail ~ 0.006 for variations of eis as small 
as Ae/s ~ 5 kJ. This gives a contribution to dfvib/dejs, 
the last term in Eq. El of at most RTAn/Aejs ^ 0.04. 
This supports our assumption that dfvib/deis = 0. 



B. Liquid Entropy 

To exploit the AG relation, we require the absolute 
value of Sc^ not just changes in Sc from one state point to 
another. To evaluate Sc{Tq) in Ea.llOlreauires S{To), the 
absolute entropy of the BKS liquid, which we calculate 
via thermodynamic integration starting from a system 
for which the entropy is known exactly. 

As our starting point we use the analytic result for 
the entropy of an ideal gas composed of two species of 
particles, each with its own mass (s^ . 



SiG = Ns^k{ln 



Nnkl In 



V f2TTmsikT 



V l2-KmokT\ ^ 



No \ 



fcln 27r0Vs^ 



(18) 



For simplicity, we continue to label the two species as "Si" 
and "O" . Note that h is Planck's constant, and that the 
Stirling approximation has been employed in the deriva- 
tion of this result, i.e. In TV! fa N\nN - N + \n{2nNy^^. 

For the purpose of thermodynamic integration, we ap- 
proximate this ideal gas with a dilute binary LJ system 
in which the stoichiometry of the species is the same as 
that of Si and O in our silica simulations, i.e. Nsi — 444 
and No = 888. 

We equilibrate the LJ system at a reference T ~ Tq — 
4000 K over a range of V from a reference V = Vq = 
8.6804 cm^/mol to V = 173610 cm^/mol. Denoting the 
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reference state point at (To,Vo) as "C", the entropy of 
the LJ Hquid S^j at C can be written, 



Slj{C) = Sig[C) + 



Ulj{C) 
T 



1 
T 



Pl-jdV. (19) 



Vo 



S'jq(C) = 1 17.236 J/mol K is calculated from Eg. [THl 
From simulations we obtain Ulj{C) — —134099 ± 
50 J/mol, the potential energy of the LJ system at C. 
Pffj = Plj — NkT/V is the excess pressure of the LJ 
system, which we evaluate for many V at Tq (Fig. |3J| . 

In order to evaluate the integral in Eq. we seek 
a function to fit to our P£j data that we can integrate 
analytically. A natural choice is the virial expansion for 
P£j, a power series in \/V . At sufficiently large the 
V dependence of P£ j will be well approximated by the 
leading term in the virial expansion, 



jjex 



b2kTN^ 



1/2 



where for our binary system. 



62 = l/9fef^^ + 4/9 6f° + 4/9&2 
and the coefficients bt^" are defined by. 



00 



(20) 



(21) 
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FIG. 5: Isotherm of |P£}| for the LJ fluid at T = 4000 K. 
The solid line is the flt to Eq. with M = 8. The dashed 
line is the curve given by Eq. 1201 The data are shown on a 
log- log plot to simplify the comparison of the data to Eq. 1201 
at large V. The cusp near V = 10 cm^/mol is due to the fact 
that P£j changes sign. 



-An 



^2t-'i>L.,(ti,v,r)/kT 



\)dr. 



(22) 



We calculate 62 numerically for the LJ system and find 
62 = —0.0337 nm^. As shown in Fig.jSJ we have simulated 
the LJ fluid to large enough V so that P£j conforms to 
Eq.|2ni To integrate P£j over (Vb, 00), we fit the data to 



pex 



b2kTN^ 

y2 



On 



M 

X! y 

n— 3 



(23) 



and use this form to evaluate the integral in Eq. ^| We 
evaluate the integral using three different fits with M — 
6, 7 and 8 in order to obtain an error estimate for the 
integral. We estimate the value of the integral to be 
-1268 ± 700 J/mol, and thus find Slj{C) = 84.028 ± 
0.175 J/mol K. 

To obtain S{C) from Slj(C) we perform a generalized 
thermodynamic integration [38(|, in which a parameter A 
is used to create a continuous path between the LJ system 
at C and the BKS system at C. To this end, we conduct 
MD simulations of a system of particles interacting via a 
pair potential $ such that, 



$(A) = A$BK5 + (1 - A)$LJ. 



(24) 



When A = 0, the system corresponds to the LJ fiuid, 
and when A = 1, the system corresponds to the BKS 
potential. For arbitrary A, the instantaneous potential 
energy is given by. 



Ux = XUbks + (1 - X)Ulj, 



(25) 



where IAbks ^lj) is the instantaneous potential energy 
of the system evaluated using only the (^bks {^lj) pair 
potential. The Helmholtz free energy difference AF — 
Fbks ~ Flj between the BKS and LJ systems at C is 
given by. 



AP(C) 




/ {UBKS-UL.j)d\. 

Jo 



(26) 



We evaluate the above integral by simulating the system 
governed by Eq. [21 at several values of A and using a 
cubic spline to interpolate between points. The integrand 
is shown in Fig. 1^1 from which we obtain AF{C) — 
— 1635990 ± 50 J via numerical integration. 

From our BKS simulations we find U{C) = -1802257± 
100 J/mol, from which we evaluate AU{C) = U{C) - 
Ulj{C). Using S{C) = Slj{C) + [AU{C) - AF{C)]/T, 
we find for the BKS liquid 5(C) = ^(ro, Va) = 75.986 ± 
0.177 J/mol K. 

To obtain values of S at Tq for V different from Vb, 
we carry out a thermodynamic integration of P along an 
isotherm of the BKS liquid, using. 



s{v,n) = s{c) + -[u{v,n)-u{c)\ 



P{V')dV'. 



(27) 
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FIG. 6: Plot of the integrand / = {Ubks - Ulj), in Eal^ 
The solid line is a cubic spline fit to the data. 
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FIG. 7: Thermodynamic properties of liquid BKS silica along 
the To = 4000 K isotherm: (a) S, (b) P, and (c) U. 



To evaluate the above integral, we find P for various V 
at To and fit the data with a cubic spHne. This spline fit 
is then used to generate data for a numerical evaluation 
of the integral. We thus have S{V,Tq), the absolute en- 
tropy of the BKS liquid at all V studied, at Tq = 4000 K 
(Fig.0. 



V (om^/mol) 

FIG. 8: U{V,T — 0) for quartz, coesite and stishovite. Solid 
lines are common tangent constructions used to find the sys- 
tem ground state energy for specific values of V (as indicated 
by dotted lines). 



C. Crystalline Ground States 

As discussed in the next section, we also find it useful 
to calculate the T = crystalline ground state energy 
U (0) of the BKS system, for comparison with the fS en- 
ergies obtained from the quenched liquid configurations. 

To this end, we study three crystalline structures of sil- 
ica important in the V range under consideration, namely 
quartz, coesite and stishovite H^ESSl- We evaluate 
the T — energy curves for these crystals as modeled by 
the ^BKS pair potential. Starting from the previously 
determined crystal structures, we optimize U{0) of the 
model system through an iterative procedure where we 
alternately minimize U (0) as a function of the particle 
coordinates in a simulation cell of fixed geometry, using 
a conjugate-gradient procedure; and then optimize the 
cell geometry with a simplex method [s^ . During the 
cell geometry optimization, we constrain V to be fixed, 
but otherwise allow the shape to change. This is done 
to remove anisotropic stress within the crystal while pre- 
serving V. Once the crystal structure has been optimized 
for a particular V, we incrementally change V and repeat 
the optimization. The results for the three crystals are 
shown in Fig.|Sl At fixed V, the thermodynamic ground 
state may be a single crystal phase, or a coexisting mix- 
ture of two crystalline phases. To obtain the ground state 
energy in the case of a mixture, we employ the "common 
tangent construction," as shown in Fig. |S| 



IV. RESULTS AND DISCUSSION 

The calculations described above yield a complete ther- 
modynamic description of the BKS model of liquid silica. 
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including the absolute free energy of the model, and the 
absolute configurational entropy, over a wide range of 
V-T conditions. Combined with dynamical data, in the 
form of the diffusion coefficient D, a number of conclu- 
sions may be drawn, as described below. Note in the fol- 
lowing that by D we mean the diffusion coefficient of the 
Si atoms, evaluated from the particle mean squared dis- 
placement in equilibrium. For comparison we also show 
some results for the diffusion coefficient of the O atoms. 
Do- In general, the qualitative results are independent 
of the choice of atomic species. 



A. Signature of fragile-to-strong crossover in the 
potential energy landscape 

Our results for D are shown in Figs. (SJa) and llOf a). 
These plots take the form of Arrhenius plots of D and 
D/T respectively, the latter quantity being preferred in 
some works as a measure of particle mobility in liquids. 
As first observed in Ref. jSj, the FSC of the BKS model 
can be seen in our simulations at low p, which correspond 
to the value of p for real silica at ambient P; isochores 
are curved on an Arrhenius plot at high T, but becomes 
straighter at the lowest T. 

We compare the behavior of D and Dq in Fig. ^] We 
find that the ratio of Dq/D varies between 1 and 3, but 
that the T dependence of Dq is qualitatively the same 
as that of D. 

Our results for e/5 and Sc are shown in Figs. ll2lll3l and 
n~n Note that our estimates for Sc have changed some- 
what, compared to the values published in Ref. [sJl, due 
to improved averaging using more data, as well as refine- 
ments in our analysis. However, the qualitative results of 
Ref. |3l| remain in agreement with those presented here. 

Consistent with the predictions made in Sect ion UTTl we 
find that for the low p isochores, where a FSC is observed, 
we also observe a point of inflection in the T dependence 
of both ejs and Sc- This inflection is what one would 
expect if the emergence of strong liquid behavior with 
decreasing T is associated with the approach of e/5 and 
Sc to a constant at low T. 

Ref. showed that for a binary Lennard-Jones mix- 
ture at low T, eis{T) — y • This T-dependence of e/5, 
consistent with a Gaussian distribution of IS energies, 
has been observed in other models 9, 24]. Furthermore, 
the binary LJ system is a relatively fragile liquid, and in 
this regard, our results for silica at high p are similar to 
those for binary LJ (inset. Fig. If 2|l . This is consistent 
with the fact that at high p the network structure of liq- 
uid silica is disrupted, giving behavior more like that of 
simpler liquids, such as the binary LJ system. We also 
find that the T-dependence of Sc becomes more like that 
of a simple, fragile liquid (Fig. If 4|l as p increases. 




1000/TS„ (mol/J) 

FIG. 9: (a) Isochores of D, the diffusion coefficient of Si 
atoms, (b) Test of the AG relation along isochores of D; the 
legend indicates p in g/cm'^ . If the data fall on a straight line, 
the AG relation is satisfied. For comparison, the inset shows 
Do, the diffusion coefficient of O atoms, along the p — 3.01 
(triangles) p — 2.31 g/cm^ isochores (stars). 



B. Implications for the Kauzmann paradox 

In Fig. El we compare the behavior of e/5 to the T = 
energy of the crystalline state of the system, as found 
from the data in Fig. |H1 It is interesting to note how 
closely e/5 approaches the crystal energy at low p, com- 
pared to the behavior at higher p. We do not expect e/5, 
the energy of a disordered configuration obtained from 
the liquid at T, to ever be less than that of the T = Q 
crystalline state of the system. The T-dependence of e/5 
at low p is consistent with behavior that would respect 
this constraint as T — > 0. At higher p, e/5 does not ap- 
proach the crystal energy as closely as it does at lower p. 
The conditions that might induce an inflection in e/5 are 
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FIG. 10: (a) Isochores of D/T. (b) Test of the AG relation 
along isochores of D/T; the legend indicates p in g/cm'^. For 
comparison, the inset shows Do along the p — 3.01 (triangles) 
p — 2.31 g/cm^ isochores (stars). The lines in the main panel 
are fits of a straight line to the data, used to obtained the 
values of A and no shown in Fig. ^] 



therefore not realized in the T range of our simulations. 

The inflection in e/5 is associated with an inflection in 
the T-dependence of Sc, also found at low-p (Fig. [T^ . 
For real systems, the third law of thermodynamics re- 
quires that the lower bound for Sc be zero. Although our 
system is purely classical, the same constraint applies, 
because the configurational entropy we calculate counts 
the number of basins explored by the liquid, which is 
necessarily one or greater. As pointed out by Kauzmann 
in 1948 |42|, the entropy as T decreases of many super- 
cooled liquids initially decreases at a sufRciently high rate 
so as to suggest that the entropy might reach zero at fi- 
nite T (the so-called "entropy catastrophe"). That this 
purely thermodynamic event seems to be preempted by 
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FIG. 11: (a) Arrhenius plot of Do (open symbols) and 
D (filled symbols) along the p — 3.01 (circles) and p — 
2.31 g/cm'' (diamonds) isochores. (b) Ratio Do/D as a func- 
tion of l/T along the isochores shown in (a). 
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FIG. 12: eis(T) along isochores. The lines are fits to each 
isochore of a fifth order polynomial in T with no linear term 
(i.e. zero slope at T = 0); these are the polynomial curves 
we use to estimate the integral term in Eq. |U] The legend 
indicates p in g/cm^. The inset shows ejs versus l/T for 
three isochores spanning the density range. The low density 
isochore shows a marked departure from the relation ejs ~ 
l/T at low T. The symbols used in the inset correspond to 
the same p as in the main panel. 
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FIG. 13: Detail of isochoric e/g behavior for three p spanning 
the range of our calculations. The thick horizontal lines show 
the value of the T = crystal energies obtained from Fig. |S| 
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FIG. 14: Sc along isochores; each panel is labeled by the 
density in g/cm"^. Each curve is obtained using Eq. |5| by 
integrating the fitted curves for eis shown in Fig. 1121 



the occurrence of a kinetic event, the glass transition, is 
the so-called "Kauzmann paradox." At low p and high 
T, we find that Sc behaves in similar fashion, decreasing 
rapidly as T decreases. An extrapolation of the observed 
high T behavior raises the possibility that Sc might reach 
zero at finite T. However, the inflection observed in the 
lower part of our observed T range establishes behav- 
ior that allows Kauzmann's "entropy catastrophe" to be 
avoided through a purely thermodynamic phenomenon. 

It is therefore tempting to speculate that our obser- 
vations may be transferable to other systems to which 
the Kauzmann paradox seems to apply. How > is 
maintained in deeply supercooled liquids can perhaps be 
understood in terms of the PES changes observed here. 
Moreover, the PES change we find in BKS silica is corre- 
lated to the fragile-to-strong dynamical crossover. Hence 
it is possible that the FSC and (the avoidance of) the 
Kauzmann paradox are fundamentally interrelated phe- 
nomena. 

However, if the above speculations are confirmed, it 
is important to note the differences between silica and 
other supercooled liquids. The picture developed above 
implies that the T range of the phenomenon by which 
silica avoids the Kauzmann paradox is above, and widely 
separated from, Tg in silica. In other supercooled liquids, 
the glass transition may occur at T above, and thus ob- 
scure, the PES changes found here. More work on these 
possibilities is clearly required. 



C. Test of the Adam-Gibbs relation 

In order to draw the conclusions given above, we must 
also test that the liquid satisfies the AG relation. If Sc 
does not control the behavior of Z?, then we will lack 



the basis required for making a connection between the 
behavior of the PES and the liquid dynamics. 

We perform this test by plotting log(£') [Fig.l^b)] and 
log(D/T) [Fig.lintb)] against log(l/r5c). The AG rela- 
tion is obeyed by data that follows a straight line on such 
a plot. We find that the isochores of D/T provide the 
best agreement with the AG relation. Note also that 
both high and low p isochores, regardless of whether Sc 
exhibits an inflection, obey the AG relation. Thus we see 
that regardless of dynamical regime (fragile or strong), 
and regardless of inflections in the T dependence of Sc or 
e/s, the liquid behaves so as to satisfy the AG relation. 
This observation reinforces the positive tests of the AG 
prediction that have been documented in other work (see 
e.g. Ref. '^,'2^121141). 

In Fig. El we present estimates of the constants A and 
/io that appear in Eq. These are obtained by fitting 
straight lines to the isochores in Fig. llOf b). omitting the 
three data points at the highest T, where deviations from 
the AG relation are expected. 

D. Entropy and Diffusion 

In the case of simulated water 20], it was found that 
maxima of Sc isotherms occur, within error, at the same 
V as the maxima in isotherms of D. We show our re- 
sults for the V dependence of Z?, Sc and S in Fig.^] At 
lower T, isotherms of Sc and D pass through a maximum 
at approximately the same V . However, as T increases, 
the correlation of these maxima fades. This may be due 
in part because at higher T, our estimates of Sc worsen, 
due to the larger role played by the anharmonic correc- 
tions. At the highest T, the trend is for isochores of Z?, 
Sc and S to become monotonic functions of V . An obser- 
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FIG. 15: Estimates of the parameters (a) fio and (b) A, that 
occur in Eq.0 These estimates are based on the straight-line 
fits to the data shown in Fig. llUf b). 
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FIG. 16: Isotherms of the V dependence of D (filled 
circles), Sc (squares) and S (triangles) at various T. 
For plotting purposes, S has been shifted down so that 
S{Vo) — Sc{Vo). The shifts in S for the various pan- 
els are (a) -66.6042 J/mol K, (b) -70.6031 J/mol K, (c) 
-77.8919 J/mol K, and (d) -84.3673 J/mol K. 
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FIG. 17: Comparison of U and the contributions to Cv- The 
top panels (a)-(c) show isochores of U spanning the studied 
density range; also shown are estimates for U of the crystalline 
state of the system in the harmonic approximation, derived 
from the T = estimate of U in Fig. |H1 and extended to 
higher T using a straight line of slope 3R/2. The bottom 
panels (d)-(f) show the contributions to Cv- The solid line is 
Cv - {3R/2){N - 1)/N- the dot-dashed line is C^''; and the 
dashed line is Ct/ . Each curve is obtained by differentiating 
the function fitted to the data for the corresponding energy. 



E. Specific heat 



In terms of the various contributions to E = ejs + 
Eanh + 5RT{N— 1) /N, we can write the constant volume 
specific heat Cv as, 



dE 
df 
d_ 
dT 



eiS + Eanh 



3RT(N - 1) 
N 



(28) 



So written, we can separately evaluate the contributions 
to Cv from e/5 and Eanh: which we denote Cy^ and C™'* 
respectively (Fig. [T7I). At all p, C^"'' exhibits a maxi- 
mum; at the lowest p the inflection in the T dependence 
of e/5 means that Cy^ also passes through a maximum. 
Together, at low p, the strength of these two contribu- 
tions becomes large enough to give a peak in the total 
value of Cv- This peak is therefore a thermal signature 
approximately demarcating the crossover from fragile to 
strong dynamical behavior. It would be interesting to 
explore if such signatures could be observed in high T 
experiments on silica, or related systems [45| . 



vation that the V dependence of entropy follows that of 
D would be consistent with recent work examining the 
relationship between structural disorder and diffusivity 
in BKS silica |4J]. 



F. Relation to polyamorpiiism 

Real amorphous solid silica displays "polyamorphism," 
the conversion under pressure of a low density form to a 
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high density form, that occurs in some ways as though 
it were a first-order phase transition. Computer sim- 
ulations of BKS siHca have provided evidence that this 
polyamorphic transition may correspond to a sub-T^ rem- 
nant of a hquid-Uquid phase transition occurring in the 
equihbrium liquid |4fil |. 

Having found that the same model, BKS silica, ex- 
hibits a thermodynamic anomaly, in the form of a Cy 
peak associated with a FSC, it is natural to ask if this 
phenomenon is related in some way to polyamorphism. 
It is difficult at present to answer this question decisively, 
since the region of the proposed liquid-liquid instability 
in BKS silica has only been approximately located, and 
seems to lie in a T range below that at which we can 
evaluate equilibrium liquid properties with current com- 
putational resources. However, several trends suggest a 
connection. 

First, we find that the T at which the peak of Cv 
occurs decreases with increasing p, and passes outside 
of our range of observation at the approximate p where 
liquid- liquid phase separation is proposed (Fig. I18|l . This 
behavior is consistent with the observed Cy peaks being 
a line (in the V-T plane) of high-T, non-singular thermo- 
dynamic anomalies that becomes singular as the critical 
region of the liquid-liquid transition is approached. 

Second, the observed line of Cy peaks naturally defines 
a "crossover zone" in the behavior of the liquid between 
a high-T, high-p region (Region I in Fig. I18() . within 
which the liquid is more fragile, the IS energies are rela- 
tively high, the tetrahedral network is disrupted, and the 
properties are in general more similar to simpler liquids; 
and a low-T, low-p region (Region II in Fig. 118(1 within 
which the liquid is becoming strong, the IS energies arc 
dropping (perhaps toward a lower limit), the tetrahedral 
network is becoming prominent, and complex thermo- 
dynamic behavior (e.g. negative expansivity) emerges. 
It is possible that these two regions of behavior, as T 
decreases, become progressively more sharply separated, 
perhaps ultimately by a first-order phase transition. 

More research, both through experiments and simu- 
lations, is required to confirm or refute such a picture. 
However, our current understanding of the BKS system 
suggests that three distinct phenomena may in fact be in- 
terrelated: (i) the FSC, (ii) polyamorphism, and (iii) the 
landscape behavior that allows the Kauzmann paradox 
to be avoided. In this regard, Sasai has recently stud- 
ied the interrelationship of the FSC and a liquid-liquid 



phase transition in a random energy model [43 ■ Also, it 
is worth considering the behavior of other systems (e.g. 
BeF2 ^^) that display one or more of the above three 
phenomena, to test if the others also appear. In partic- 
ular, for any system that becomes strong at low T via 
a FSC, it may be that polyamorphism can be observed 
under nearby thermodynamic conditions. This may be a 
useful clue for identifying polyamorphic materials. 
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FIG. 18: Location of the line of Cv maxima (asterisks) in 
the V-T plane. Also shown are points on the "temperature of 
maximum density" (TMD) line (triangles), at which the iso- 
baric expansivity changes sign; and the location of maxima 
in the contribution of e/s to Cv (squares). The diamond in- 
dicates where evidence for liquid-liquid phase separation was 
found in Ref. li^ . The regions labeled I and II are referred 
to in the text. 
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